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We propose new classes of random matrix ensembles whose statistical properties are intermediate 
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based on integrable A'^-body classical systems with a random distribution of momenta and coordi- 
nates of the particles. The Lax matrices of these systems yield random matrix ensembles whose 
joint distribution of eigenvalues can be calculated analytically thanks to integrability of the under- 
lying system. Formulas for spacing distributions and level compressibility are obtained for various 
instances of such ensembles. 
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I. INTRODUCTION 



The theory of random matrices, introduced by Wigner in the 1950s, has proved to be a very useful tool in many 
fields of physics, from localisation theory to quantum transport (see e.g. [1] and references therein). In quantum 
chaos, a well accepted conjecture states that Wigner-Dyson random matrix ensembles describe statistical properties 
of spectra of quantum systems whose classical counterpart is chaotic 0] , while statistics of integrable systems is best 
described by Poisson statistics of independent random variables [3] . The corresponding wave functions are extended in 
the chaotic case and localised in the integrable case. The choice of the random matrix ensemble suited to describe the 
statistical behaviour of a system depends on the symmetries of that system. In the usual setting [3] , standard random 
matrix ensembles consist of matrices M with independent Gaussian random elements whose measure is invariant over 
conjugation 

M — >U'^MU, (1) 

where U is an arbitrary matrix belonging to one of the three following groups of matrices: unitary, orthogonal, or 
symplectic. The unitary group defines the Gaussian Unitary Ensemble (GUE), which is supposed to describe statistical 
properties of energy levels of chaotic systems without time-reversal invariance. The orthogonal group corresponds 
to Gaussian Orthogonal Ensemble (GOE), used for time- reversal invariant chaotic systems. The symplectic group 
gives rise to Gaussian Symplectic Ensemble (GSE), applicable to time-reversal chaotic systems with half- integer spin 
without rotational symmetry. 

Though many extensions and generalisations of random matrices have been proposed in order to best describe 
various models [5|, the existence of a large invariance group as in ([T]) remains their characteristic feature. Without such 
invariance it is very difficult to connect analytically simple properties of matrix elements with complex properties of 
matrix eigenvalues. For all random matrix ensembles with invariance group it is possible to integrate over unnecessary 
variables in order to get explicitly the joint distribution of eigenvalues under the form 

P(Ai, . . . , Ajv) ^ n I^J- - ^'l^e' (2) 

with V{x) a system-dependent potential and /3 a parameter. For the Gaussian ensembles the potential is quadratic and 
the parameter /3 is equal to 1 for GOE, 2 for GUE, and 4 for GSE. All correlation functions for invariant ensembles can 
be calculated analytically [^]. However the resulting formulas are cumbersome. For the nearest-neighbour distribution 
P(s), instead of the exact expression one often uses a simple surmise proposed by Wigner. This surmise has correct 
functional dependence at small and large argument and takes the form 

P{s) = as^e-^"' (3) 
with constants a and h determined from the normalisation conditions 

I P{s)ds = / sP{s)As = 1 . (4) 
Jo Jo 

The Wigner-type surmise for the probability P{n, s) that between two eigenvalues separated by s there exist exactly 
n — 1 other levels (with P{1, s) = P{s)) is [q| 

P{n, s) = ans'^-e-^-"^ , d„ = n - 1 + ]^n{n -t- 1)/3 (5) 
and On, hn are fixed by the normalisations 

I F(n,s)ds = l, / sF(n,s)ds = n. (6) 
Jo Jo 

While for chaotic systems it is possible to argue that eigenstates may statistically be invariant under rotations, 
this is not the case for more general models. In order to describe statistical properties of such systems one has to 
consider non-invariant ensembles of random matrices. One of the most investigated examples is the three-dimensional 
Anderson model [tI], with on-site disorder and nearest-neighbour coupling. Depending on the strength of the disorder, 
it can display metallic behaviour well described by standard random matrix ensembles, or insulator behaviour with 
Poisson-like spectrum. However, at the metal-insulator transition, spectral statistics are of an intermediate type and 
are not described by invariant ensembles [8| . Similar behaviours have been observed in pseudo- integrable billiards [QJ , 
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quantum maps corresponding to difFractive classical maps [10|, or quantum Hall transitions Models have been 
proposed to describe such intermediate statistics [12j |. and random matrix ensembles which possess similar features 
have been constructed, e.g. power-law random banded matrix ensembles [isl flU. 

The main purpose of this paper is to construct random matrix ensembles which are not invariant over rotations of 
eigenstates, but whose joint distributions of eigenvalues can nevertheless be calculated analytically. A short version 
of the paper has been published in [l5|. All of these ensembles have intermediate statistics, and for certain of them 
spectral correlation functions, e.g. the nearest-neighbour distribution, are obtained explicitly. Eigenfunctions of these 
ensembles are neither localised (as for integrable systems) nor extended (as for chaotic models) but have fractal 
properties [l6| . 

Random matrices of the proposed critical ensembles are constructed from the Lax matrices of classical integrable 
models. These models are systems of N classical particles labelled by an index i, 1 < i < N in a, one-dimensional 
space. Each particle i is characterised by its position in space qi and its momentum pi. The dynamics of the particles 
is entirely described by the Hamiltonian i7(p, q), where p — {pi, . . . ,pjq) and q = (gi, . . . , ^at). The characteristic 
property of these models is the existence of a pair oi N x N matrices L and M, called the Lax pair of the system 
such that the equations of motion (the Hamilton equations, derived from the system Hamiltonian) are equivalent to 

dL 

— ^ M L- LM. (7) 

The Lax matrix L is a matrix depending on momenta p and coordinates q. We propose to consider these Lax matrices 
as random matrices with a certain 'natural' measure of random variables pj and qj 

dL = P(p,q)d^pd^q . (8) 

The explicit form of this measure depends on the system and will be discussed below. We do not impose any dynamics 
on variables p and q. The only information we use from the integrability of the underlying classical system is the 
existence and explicit form of action-angle variables /^(p , q) and (/)q(p , q). In particular, it is well known that the 
transformation from momenta and coordinates to action-angle variables is canonical, so that 

]Jdp,d9j =[|d/„d(^o . (9) 
j « 

Direct proof that the transformation is canonical is difScult in general, and implicit methods have been used to establish 
it for specific systems 18]-[2flj. In the models we consider here, action variables turn out to be the eigenvalues Aq of 
the Lax matrix, or a simple function of them. The canonical change of variables from momenta and coordinates to 
action-angle variables in ([5]) leads to a formal relation 

dL = 7'(A,(/))d^Ad^0 , (10) 

where V{\ , ) = P(p(A , ) , q(A ,</>)). The exact joint distribution of eigenvalues is then obtained by integration 
over angle variables, which can easily be performed in all cases considered, and yields 

P{\) = j 7'(A,0)d^0. (11) 

This scheme is general and can be adapted to many different models. 

In this paper we consider in detail four typical models of A^-particle classical integrable systems. The three first, 
labelled CM^, CM/j, and CMj, correspond to the rational, hyperbolic, and trigonometric Calogero-Moser models 
[2l|,[23]. The fourth model, labelled RS, is a trigonometric variant of the Ruijsenaars-Schneider model (23j . 

The Calogero-Moser models are defined by the Hamiltonian 



j<k 



where u(^) is a potential depending on the distance between particles and is a constant |24 |. For the models 
considered here it has the form u(^) = x^iO^ where 



<0 



sinh(/i^/2) 
m/2 
. sin(/i^/2) 



model CMr 

model CMh . (13) 
model CMt 



4 



The Hamiltonian of our fourth model, the trigonometric Ruijsenaars-Schneider model, is [2S 



N 



iI(p,q)=^cos(p,)n 1 

j=i k^j \ 



^5/2] 



1/2 



model RS 



(14) 



The plan of the paper is the following. Sections |TI1 IIII[ and IIVI are devoted to the construction of critical ensembles 
related respectively with the rational, hyperbolic, and trigonometric Calogero-Moser models. In each of these sections 
we briefly present the construction of the action-angle variables and choose a 'natural' measure of random momenta 
and coordinates which allows an easy change of variables as in (jlOp. We then give explicit formulas for the joint 
distribution of eigenvalues for the resulting critical ensembles of Lax matrices. In section |V] this scheme is applied 
to the Ruijsenaars-Schneider model. For this model the joint distribution of eigenvalues takes a form which makes 
it suitable for the application of the transfer operator formalism. This approach is detailed in section IVI[ and in 
section IVlIl it is applied to the analytic calculation of nearest-neighbour distributions for the RS model. The spectral 
compressibility for this model is obtained in section IVIIIl 

For clarity we state below the principal results for the four models considered in this paper. 



CM,, ensemble 

The CMr ensemble is defined as the ensemble oi N x N Hermitian matrices of the form 

^kr = PrOkr + Iff , 

qk - qr 

with g a real constant. Positions q and momenta p are random variables distributed according to the density 



(15) 



P(p,q) ~ exp 



(16) 



with A and B arbitrary positive constants. The joint distribution of eigenvalues for this ensemble is then given by 



P(A) ~ exp 



(17) 



A characteristic property of this ensemble is the exponentially strong level repulsion: the nearest-neighbour spacing 
distribution P(s) is characterised by 



In P(s) 



0(1). 



(18) 



We propose the following Wigner-type surmise for the next-to-nearest-neighbour spacing distributions P{n,s), de- 
pending on four parameters: 



(19) 



P(n, s) — as exp — 2 ^ 



It contains two fitting constants depending on n. The other two are fixed by the normalisation ([5]). 



CMh ensemble 

The CM/i ensemble is defined as the ensemble of x iV Hermitian matrices of the form 

Lkr = PrOkr + Iff „ . , r , rjTT (20) 

2sinh [^{qk - qrj/2\ 
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with g and fi real constants, and q and p distributed according to the density 



P(p,q) - exp 



'^^?''^'^''S4sinh^[M<Z,-'Z.)/2] 



— B cosh /igj 



(21) 



The exact joint distribution for this model is 



P(A)^exp -A^A^ mifo Un |l 



(22) 



where K{)[x) is the modified Bessel function of the second kind. The nearest- neighbour spacing distribution has an 
exponential asymptotic similar to (|18p but with 1/s leading term instead of namely 



The Wigner-type surmise for CM^ is 



In P(s) - — + 0(lns). 

s->0 S 



P{n, s) = as'* exp ^ cs^ 



(23) 



(24) 



CMt ensemble 

Matrices from this ensemble correspond to a situation where fi in Eq. (1201) is allowed to take pure imaginary values. 
They are of the form 



Lkr = PrSkr + 157 



/i(l - Skr) 



2 sin [^{qk - qr)/2\ 

with g and fi real constants, and q and p distributed according to the density 



P(p,q) - exp 



, ^2 



(25) 



(26) 



with the restrictions that all qj are between and 27r//i. The exact joint distribution of eigenvalues for this ensemble 
is 



P(A)~exp|^-AE^') X(A) 



(27) 



where the function x(A) is equal to 1 if Ai < A2 < . . . < Xn and A^+i — Aq > fig for all a. The nearest-neighbour 
spacing distribution is given by a shifted Poisson distribution of the form 



0, 



< s < 6 



Pis) 



1 f 

exp — , s > b 



1-6 



1-5 



(28) 



with b some fitting constant. 



RS ensemble 



The RS ensemble is defined as the ensemble of N x N matrices of the form [23 

sin [/igcr/2] 



sin U(qfc - qr + gcr)/2] 



.yi/2ei-P,V2 



(29) 
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with 

_ TT sin [fijqk - qj - ffcr)/2] ~ _ sin - + ffg)/2] 

^'^-n sin[Mfe-<Z,)/2] ' ^'=-1,1 sin - ,,)/2] ' ^'"^ 

Let ft be the set of q such that for all k the sign of both Vk and Wk is the same as the sign of sm{N figa / 2) / sm{iig<T /2) . 
The matrix L is unitary if and only if q e 17. The variables q and p are chosen to be distributed according to the 
uniform density in the region where L is unitary. That is, we choose momentum variables pj independent and 
uniformly distributed between and 27r/cr and coordinate variables q uniformly distributed over fl. In this case 
eigenvalues of the Lax matrices (j29p are also uniformly distributed over Q. Choosing ~ 27r/iV, cr = 1 and g ~ a, we 
compute correlation functions of eigenvalues of matrix (|29p for fixed a and N — >■ oo. The results strongly depend on 
the integer part of a. For < a < 1 the nearest-neighbour spacing distribution is similar to (|28p with constant b now 
equal to a. For 1 < a < 2 the nearest-neighbour distribution takes the form 



^^sinh^(ps) when 1 < 5 < 4/3 

81 „2 
64* 

sin (ps) when 4/3 < 5 < 2 



P(s) = { |is2 ^Yien g = 4/3 . (31) 



Constants A and p are determined from the normalisation conditions Q . Other correlation functions are also obtained 
in section IVlIl 

Numerical implementation 

The results presented above are quite robust with respect to alterations in the distribution of q and p. In all models 
considered we chose (as explained above) a distribution of coordinates such that the qj are confined to a finite interval 
while having a strong repulsion between each other. As may be expected physically (though we do not have a rigorous 
proof for this), numerical evidence shows that, if we keep these two characteristic features, spectral properties for 
N ^ 00 depend only weakly on the precise choice for the distribution of q and p. 

From these considerations it is thus natural to use, rather than the exact complicated distribution of q, the picket- 
fence configuration when all coordinates are just fixed and equally spaced. As all definitions of our ensembles involve 
only differences qj — qu multiplied by a parameter (/i or g, depending on the model), we can without loss of generality 
choose to take qj — j , j — 1, . . . N . 

For numerical implementation we chose qj — j, and pj as independent Gaussian variables with zero mean and with 
variance equal 1 (CM ensembles) or independent variables uniformly distributed between and 2tt (RS ensemble). 
For concreteness, we fixed n = A-k/N for CM;i and CMf and ^ = 2'k/N for RS. For such a choice, the N x N Lax 
matrices take the form 

1 — 5kr 

Lkr = PkOkr + i.g-^j — model CMr 
27r(l - 5kr) 
TVsmh 27r(fc-r)/ArJ 

2tt{i -5kr) ■ ^^^> 

Lkr=Pk Skr + 15 . r„ ^TTTT , modcl CMt 

N sm 1 27r(fc ~ rj/N\ 
1 — 

Lkr — e'^*" — -, . , model RS 

For CMt matrices with even iV, to avoid the singularity we changed — ^ iV + 1 in the above formula. As the figures 
in the next sections show, despite this particular choice for the distribution of q and p, the agreement between the 
computed spectral statistics and analytical formulas is remarkable. 

II. RATIONAL CALOGERO-MOSER MODEL 

The first model we consider is the rational Calogero-Moser model CM^. [23] , characterised by the Lax matrix 

Lkr = PrSkr + 15 — ■ (33) 

qk - qr 
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It depends on a real constant g and on a set of 2A'' random variables pk and whose distribution will be specified 
later on. We are interested in eigenvalues Aq, and eigenfunctions Ufc(Q;) of this matrix (here and below we will use the 
Greek letters to label eigenvalues and corresponding eigenfunctions) 



JV 



r=l 

To construct angle-action variables let us define the new quantities 

Qa/s = ^ul{a)qkUk{j3). 



(34) 



(35) 



From Eq. (|33|) one gets 

Lkr^r — IkLkr — " 15(1 — Skr) ■ 

Multiplying both sides by Uk{a)ur{(3) and summing over k and r one gets 

QapiK - A/3) -ig(e*e^ - Sa/i), 

where 

Cq = ^Ufe(a) . 



(36) 
(37) 
(38) 



For a = P, Eq. ([37)) implies that |eQ,|^ = 1, and one can choose the overall phase of the eigenvector Mfc(a) in such a 
way that e« = 1. Let (j)^ be new variables defined by 



Then from Eq. d^T]) we have 



Qal3 — (ka^afi — \g 



/3 



Aq — 



(39) 



(40) 



The matrix Q can be seen as the dual matrix of L, with ^q, playing the role of momenta and Aq the role of positions. 
In [l^ it was proved that there is a canonical transformation from position and momentum variables {qk,Pk) to action 
and angle variables (Aq., (/)q). Showing that the transformation is canonical is a rather technical mathematical result. 
However one can easily check that the new variables Aq, and 4>a verify Hamilton- Jacobi equations (see Appendix \K\ . 

We now consider an ensemble of Hermitian matrices of the form (I33p with random variables pk and qk drawn 
according to the measure 



P{L)dL ^ JV exp 



dpk dqk 



(41) 



J k 



where A and B are given constants and J\f a normalisation factor. The first term in Eq. (|4ip is the analog of the usual 
Gaussian weight of RMT; the second term is a quadratic confinement potential. Since the action-angle transformation 
is canonical one has 



Y[ dpk dqk = Y[ dXadct)a ■ 



(42) 



From Eq. (j35p , using orthogonality of eigenvectors one gets 

TrQ^ = E«: 



(43) 



Using these relations one can rewrite the distribution (j4ip in action-angle variables Aq and (f)a as 

1 



P{L)dL = Af exp 



(Aq - A/3)2 



dXa d(j)a 



(44) 
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Integration over the 0q gives a constant. We thus obtain the joint distribution of eigenvahies for the ensemble of 
random matrices L with the measure (j41l) as 



P(Ai, . . . , Ajv) ~ exp 



(45) 



Note that, similarly as in the standard RMT case, this joint eigenvalue distribution can be interpreted via the Coulomb 
gas model as the partition function of an ensemble of particles on a line, here with inverse square repulsion. After 
rescaling Xk = \k{Bg^ / A)~'^/^ , equilibria positions of the particles at positions Aq, are given by 

^fe = ^E (, „^,)3 .l<fc<^- (46) 

Such a relation characterises the zeros of Hermite polynomials of degree N (see also Eq. (10.3) of [24]). It is known 
from RMT Q that the distribution of eigenvalues of Gaussian random ensembles has a similar property, which implies 
that the asymptotic density of eigenvalues is given by Wigner's semi-circle law. 

An immediate consequence of the distribution (j45p is the unusual very strong level repulsion at small distances. For 
all standard random matrix ensembles the nearest-neighbour distribution P{s) behaves as at small s. By contrast, 
in our case it follows from (|45|) that 

P{s) ^ ae-''^\ (47) 

As the potential between eigenvalues decreases as the inverse square of the distance between them, the probability 
of having a gap of size s for large s is exponentially small. We could not calculate exactly correlation functions 
for the distribution (j45l) . However, combining the two asymptotics above, we build a Wigner-type surmise for the 
nearest-neighbour spacing distribution of the form 

P{s) = ac-''/'^-'"' , (48) 

where 6 is a fitting constant, and constants a and c are determined from the normalisation conditions (|4]). For the nth 
nearest-neighbour spacing distributions P{n, s), with n >2, we conjecture, by analogy with the Wigner surmise ([5]) 
for standard random matrices, the form 

P(n,s) = asV/'^'-"" (49) 

with two fitting constants h and d. 

To assess this conjecture we compare the analytical expressions (|48|) -(|49 | with numerical results, with the choice 
of parameters detailed in section HI Results displayed in Fig. [l]show that the agreement is remarkable. 

III. HYPERBOLIC CALOGERO-MOSER MODEL 

The Lax matrix for the hyperbolic Calogero-Moser model CM/i reads [l^] 

Lkr = PrSkr + ^9(1 - 4r) „ ■ i^ TTT^ ■ (50) 

2smh[fi[qk - 9r)/2) 

Let us define two matrices Q and R by 

Qcp = E 4(«)e^'^'=Ufc(/3) , Rc.p=Yl <(a)e^^*«fc(/3) , (51) 

k k 

and two vectors 

= ^ Wfc(a)e'''?''/2 , /a = ^ Wfc(a)e-''«'=/2 . (52) 

k k 

From ()50p one can get the two equivalent equations 

e^^'i^Lkr - Lkre^^"- = i.gM(l - 4r)e'^(«'=+«'-)/2 (53) 
e-'^'^-Lfc, ~ Lfc,e-^«'= = 15^(1 - <5fc,)e-'^(«''+«'^)/2 . (54) 
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FIG. 1. (Color online). Nearest-neighbour distributions P(n,s) for the random matrices of model CMr with g — 0.05 (top 
left), 0.25 (top right), 0.5 (bottom left), and 1 (bottom right), averaged over the central quarter of the unfolded spectrum for 
8000 realisations of matrices of size A'^ = 512. Solid lines are numerical results, dashed lines indicate the fits (|48|) for and (|49|) 
for P{n, s) with (in each panel from left to right) P(s) = -P(l, s) (red), P{2, s) (green) and P(3, s) (blue). 



Multiplying both sides by Uk{a)*Ur{l3) and summing over all k and r one gets 



QapiK - A/j) = -ig^(e*e^ - Qa/s) 

RapiK - A^) = iffAil/a//? - Rap) , 



(55) 
(56) 



which implies that matrices Q and R take the form 



A/3 - Aq + igfi 



Ra0 - la 



Xp - \a - ^gfj^ 



fp 



(57) 



By their definition (j5ip . matrices Q and R are inverse of each other, so that QajRjp = ^ap for all a, (3. For a 
this condition implies that 



2 2 * r \ ^ 

.9 M e„/a 2^ 



(A^ - Aq + i5/i)2 



1 



(in particular, it follows that all and are nonzero). For a ^ (3 one obtains 



E 



(A^ - Aa + i.gM)(A^ - A/3 - ig^) 







Using the identity 



(A^ - Aq + igfJ,){Xj - \p- iffAi) 
valid for a 7^ /?, one concludes that 

E 



A/3 - A^ + ig^ A/3 - Ac + 15/^ 



^7/7 



Aq, — A.-) 



A^ - Aq + ig^ 



(58) 



(59) 



(60) 



(61) 
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where c is a certain constant independent on a. According to this equation the quantities = e^/*/c obey a system 
of linear equations of the form 



E 



^7 Va 



= 1 , 



(62) 



with = and Ua = — i^M- This equation coincides with Eq. (|B1|) in Appendix [Bl From (|B2p it follows that 



where 



Aa — 



while (|B4p implies that 



(63) 



(64) 



(65) 



It readily follows from the definition (j52p of Cq and fa that ea/a = A, thus the value of c is fixed by ig^ic = 1. 
Equation ((58)) is then fulfilled as a direct consequence of (|B3p . Finally we have 

(66) 



Let be new variables defined from diagonal elements of matrix Q by 
Then from ^ and dM]) it follows that 



(67) 



(68) 



In the definition (j52p of Cq, it is convenient to choose the overall phase of the eigenvector Uk (a) in such a way that Cc 
be real. As Qaa — leap one has 



e„ = |14|i/2eM0»/2_ 

Using Eq. (|57p . the matrix Q can now be expressed in terms of the new variables Aq and 0q, as 



A/3 - Ac + ig^ 



(69) 



(70) 



As in the case of model CM^, the matrix Q can be seen as the dual matrix of L. Indeed, Q coincides with the Lax 
matrix of the rational Ruijsenaars- Schneider model with coordinates Aq and momenta (f>a [3- In [isj it has been 
proved that the transformation from position and momentum variables {qk,Pk) to action-angle variables (Aq,(/)q) is 
canonical. Again one can check that the new variables Aq and (pa verify Hamilton- Jacobi equations. 

We now consider an ensemble of Hermitian matrices of the form (I50|) with random variables pk and qk drawn 
according to the measure 



P{L)dL = J\f exp 



-ATxL^ cosh 



dpk dqk 



(71) 



J k 



As in the case of model CM^, Eq. (|71l) contains a standard RMT Gaussian weight and a confinement potential which 
can be rewritten as TrQ -I- Tri?. Using (p7|) . (1551) . ^^.d the fact that the transformation is canonical, we get the 
distribution in terms of the new variables Aq, and 6a as 



P{L)dL^Mexp 



- A Xl -BY^\Va\ cosh ficPa 



(72) 
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The joint distribution of eigenvalues is then obtained by integrating over the angle variables, using 

f°° 1 

/ exp [-B\Va\ cosh ^l<j)a]d<j)a = -Ko{B\Va\) 



(73) 



where Kq is the modified Bessel function of the second kind. This yields the joint distribution of eigenvalues for model 
CM^ as 



\ a J a y I3=ia 



Aa — A^ 



(74) 



This expression is exact but difficult to handle. In order to find a Wigner-type surmise for the nearest-neighbour 
distributions we consider the limiting behaviour P(A) when two nearby eigenvalues Ai and A2 get close to each other. 
Setting s Ai — A2 we see that the factor exp(— in the case of model CM^ is replaced by a factor 



Kn B\ 1 



s exp 



We therefore expect the nearest-neighbour spacing distribution to behave as 

P{n, s) = as"^ exp(— 6/s — cs) . 



(75) 



(76) 



In Fig. [5] we show the results of numerical computations of the nearest-neighbour spacing distributions for matrices of 
the form (|50p with the choice of parameters and variables detailed in section ID The surmise (|76p perfectly reproduces 
numerical results. 




FIG. 2. (Color online). Nearest-neighbour spacing distributions P{n, s) for the random matrices of model CM^ with /j, — An /N 
and g — 0.05 (top left), 0.25 (top right), 0.5 (bottom left), and 1 (bottom right), averaged over the central quarter of the 
spectrum for 32000 realisations of matrices of size TV = 256. Solid lines are numerical results, dashed lines indicate the fit (|76p 
for P{n, s) with (in each panel from left to right) P(s) = P(l, s) (red), P(2, s) (green) and P(3, s) (blue). 



In order to assess better the validity of the exponentially strong level repulsion for models CM^ and CM/j, we 
compare in Fig. |3] the beginning of the distributions P(s) for these models. Clearly the 1/s^ repulsion for CM^ and 
the 1/s repulsion for CM;i fit numerical curves very well. However, the precision of our numerical results does not 
permit to confirm or reject the presence of the logarithmic term dins in P(s) for CM^ model. 
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FIG. 3. (Color online). Nearest-neighbour spacing distributions P(s) for the random matrices CM^ {g — 0.5, black circles) 
and g — 1, red squares) and CM^ (g — .5, green triangles up, and g = 1, blue triangles down), with /i — Att/N. Symbols are 
numerical results, solid lines indicate the fit + cs — In a (model CMr, Eq. (|48[) ) and 6/s + cs — In a — dins (model CMh, 
Eq. (|76p ). Logarithm is natural. 



IV. TRIGONOMETRIC CALOGERO-MOSER MODEL 



For the trigonometric Calogero-Moser model CMj the Lax matrix is [24 



Lkr = PrSkr + 1.9(1 - 4r) 



2sm{p{qk ~ qr)l'2) 



(77) 



The only difference with the previous model, (|50p . is the sin function which replaces the sinh. The matrix ()77p can be 
obtained from (j50|) by the substitution /i — > i/i. However the fact that positions of the particles are now defined on a 
circle (because of the sin function) makes the resulting spectral statistics entirely different from the previous models 
CM^ and CMh. 

To construct action-angle variable we introduce, as in the previous section, two matrices 



k k 

and two vectors 

ea=5]iifc(a)e'^''-/' 

fc k 

Following the same steps as above with fi replaced by i/i, one gets 

.9Ai 



/a = ^Wfe(a)e-'''''- 



/2 



]al3 — fa- 



Again, using the fact that Q is the inverse of R we obtain that 



Rap — 



A/3 - Aq -f gfi 



E 



A'v — Aq — 



= ci , 



E 



\.f,\' 



A7 — Aq 



5M 



= C2 



(78) 



(79) 



(80) 



(81) 
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with certain constants ci and C2 independent on a. Repeating the same arguments as in the previous section and 
using resuhs of Appendix IB] one concludes that C2 — — Ci = l/{fig) and 

|ea|' = K., = (82) 

with 

The new variables 4>a are defined as above from the diagonal elements of matrix Q as follows 

Qaa = V^/^W^/^'^'^- . (84) 

In the definition (17^ of one has a freedom to choose the overall phase of the eigenvector (a) . Since from (jSU]) 
one must have Qaa — fa^a, one can choose phases, for example, as follows 

e„ = yV2eiM*«/2 ^ ^ ^V2e-i^*=/2 . (85) 

Then the matrix Q can be expressed in terms of new variables Aq and (jja as 

Q^p = e'^*°/2w^„i/2- yi/2giM0,/2 _ (gg) 

The inverse matrix i? plays a symmetric role, as it is obtained from Q by exchanging /i to — Again, there is a 
canonical transformation from position and momentum variables {qk,Pk) to action and angle variables {Xa,4>a) |l8| . 
An important consequence of ((82)) is that for all a we should have 

y„ > 0, > . (87) 

These inequalities impose non-trivial restrictions on eigenvalues A^, as we will see now. We label eigenvalues so that 
Ai < A2 < ... < Aat , and we consider the function 

Mx) = Er^ + -- (88) 

Xy-x fig 

It has N poles at a: = A^, and according to Eqs. ([5T|) and (1551) it has N zeros at x = A^ + fig. If all numerators Va 
are positive then the derivative of h is positive, and it is easy to check from the graph of the function that between 
two consecutive poles there is one and only one zero. Suppose fig > 0. Then for x — t- —00 the function h{x) has a 
strictly positive limit. The lowest zero Ai + fig must thus lie in the interval ]Ai, A2[. More generally one must have 
Aa + g]Aq, Aq-)-i[ for 1 < a < iV — 1, while the largest zero A^v + fig lies in the interval JAat, cx)[. Thus eigenvalues 
fulfill the inequalities 

Aq+1 - Xa> gfi . (89) 

Conversely, if these inequalities are fulfilled then trivially all Va are positive. Therefore, are the necessary and 
sufficient conditions for the positivity of all Va- In particular, eigenvalues of the Lax matrix (j77p obey inequalities 
([89)1 for any choice of the q. These results adapt straightforwardly to the case where fig is negative. 

Since in (j77p the qk only appear as an argument in the sin function, there is no need to choose a confining potential 
for the particle distribution as in models CM^ and GMh ■ We consider the probability distribution of p and q in the 
form 



P(p , q ) - exp 



(90) 



with the restrictions that all qj are between and 27r//i. Since the change of variables from p and q to Aq, and (j)a is 
canonical and the restrictions (|89p do not depend on phase variables the joint distribution of eigenvalues is 



P(Ai,...,Ajv)^expl -AVA^ W(Ai,...,AAr) , (91) 
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where the function x(A) is equal to 1 if (l89l) is fulfilled for all a, and otherwise. 

It turns out that model CM; is very similar to a fourth model, the Ruijsenaars- Schneider model, that we will consider 
in the next section. Therefore we postpone analytical calculations of the nearest-neighbour spacing distributions to 
section IVTl The nearest- neighbour spacing distributions P(n,s) are shifted Poisson distributions of the form 



0, 



P{n,s) 



{s - nby 



(n - 1)!(1 - 6)'^ 



,-{s-nb)/{l-b) 



< s <nb 
s > nb 



(92) 



with some numerical constant b. In Fig. 2] we show the results of numerical computations for matrices of the form 
([77)) with the choice of parameters and variables detailed in section H] 
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FIG. 4. (Color online). Same as Fig. |2] for model CMt with P{n,s) the shifted Poisson distribution (p2)) . 



V. RUIJSENAARS-SCHNEIDER MODEL 



The Calogero-Moser models considered in the previous sections are such that there exists a matrix Q which is, in a 
certain sense, dual to the Lax matrix L. Namely, the canonical transformation from variables {pk,Qk) to action-angle 
variables (A^, 0q) is such that the action variables Aq, are eigenvalues of L and angle variables (jfa are related to Qaa 
in a simple way. In fact, the matrix Q is also a Lax matrix, corresponding to a possibly different Hamiltonian, and L 
plays the role of a matrix dual to Q for the inverse of the canonical transformation [3] • The rational Calogero-Moser 
system is self-dual since matrices L and Q, given by and PO]) . are equal up to labelling of the variables. 

The model we consider in this section is related to the above models in that its Lax matrix L and the dual matrix 
Q arc a kind of generalisation of those of model CMt ([55]) . The treatment of this model closely follows the previous 
section. 

The Lax matrix for Ruijsenaars-Schneider model is the N x N unitary matrix given by [20| 

Lur ^ d^^'-l^wl'' . . f " ^ ^ . v;i/2ei<^p./2 . (93) 

sm Yn{qk - gr + 5f^j/2J 

with 

t7 TT ^™ [^(* - qj - 9cr)/2] ~ yj sin [ii{qk - q, + 9(j)/2] 

'''^'= = 11 \7^-\ — ' ^^fe = 1 1 ^-r-, rj77\ — ■ (94) 

J^k sm [/x(gfc - gj)/2j sm [^(q^ - gj-)/2j 
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This Lax matrix is related with the Hamiltonian ([14]) by 

i/(p,q) = iTr(L + Lt) . (95) 

It is convenient to introduce the vectors 

gfc = V;l/2gi^Pfc/2giM9fc/2 ^ ^ ^l/2gi^p,/2g-i,*9fc/2 ^ (-gg^, 

SO that matrix L can be rewritten as 

sm [^(gj; - gr + ffo-)/2J 

As in the previous section, the condition that the matrix (|93p is unitary imposes certain restrictions on the coordinates 
q, which we will discuss later. Assuming that the Lax matrix is unitary, we choose to denote its eigenvalues by e"^'*'" . 
The dual matrices for the Ruijsenaars-Schneider model are defined by [20] 

Qcp = J2 KW^'^'MP) , Rcp ^Q}>a=J2 <{a)e-''""'MP) , (98) 

k k 

and the vectors Cq. and fa by 

ea = ^Ufc(a)efc , fa = ^Ufc(a)/fc . (99) 
k k 

From Eq. ([^ one has 

giM(g.+9.)/2^^^ Sin [ii{qk - + ga)/2\ = /,! sin [^lga/2] . (100) 
Multiplying this expression by ul.{a)ur{p) and summing both sides over k and r leads to 

^Qcp (e-^^+-^s/2 - e-^°--^s/2j ^ ^* [A..ga/2]e^ , (101) 
which yields the analogue of ([50]) . 

Q.. = /:e---/- . r.r^T^!^ V2f -'"-^^-^ • 
sm [ct(A;3 - Aa +gM)/2J 

Let us rewrite 

e^''^ — 1 

where we have set r = cr/ig. From its definition (jOSp it is clear that Q has to be an unitary matrix, i.e. Qa-yQp-^ = 
i5q,0. Selecting terms with a = /3 and a ^ 13 yields the two equations 



iH^i/"i^E u(..„-':7Li.A.i2 =i^ (104) 



7 

|2 



[tlJ ^0 (105) 



^gi((TAtv-r) _ gi(TA^^^g-i((TA^-T) _ icrA^ ^ 

where we have set p = e"'"^ — 1. Using the identity 



(giCo-Ac— r) _ gicrA^ ^^g— i((TA^— r) _ g-icrA^, 
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it follows from (|105p that there exists a constant c such that 

7 

Equations (|B1[) -(|B2 |) allow to obtain 

leap = -cpVae 
while from Eq. (I104p one obtains, using Eq. (|B3|) 



-ir(Ar-l)/2 



|/„|2 = _lp^„e-(A'-i)/2 



with new vectors and defined by 

^ _ TT sin [g(Aa - A/j + gA^)/2] 
sin[a(A„~A^)/2] 



sin [g(Aa - Aff - g[i)li\ 
sin [ct(Aq - A/j)/2] 



Using (jBSp we get 



5:|e„p = c(l-e-'^^). 



(107) 



(108) 



(109) 



(110) 



(111) 



There is some overall freedom in the definition of vectors and in Eq. (1961) . as one could multiply Cq by some 
constant factor and divide by the same factor. This in turn entails the same freedom for vectors and in (j99p . 
If one chooses for instance |ejp = t then from unitarity of the transformation in (|99p one has l^aP = ^, which 



fixes the value 



and thus 



1 -e- 



tN 



sin(7VT/2) 



isin(r/2) 



By definition t = |ejp is positive. A convenient choice is to take 



sin(iVT/2) 



sin(r/2) 



(112) 



(113) 



(114) 



Since \ea^ and are non-negative one concludes that Vq and Wa have the same sign as sin(A^T/2)/ sin(T/2) for 
all a. The choice (11141) implies that 



led' = \Vo\ , l/aP = IW^aj . 

The new action variables are the Aq, and angle variables (f>a are defined by 



(115) 



(116) 



Equation (|102p implies that Qqq, = fa^a^ laX^ _ jj^ view of (|115p and (|116p one can choose the phases of the 
eigenvectors Uk{a) such that and /q. defined by (j99p can be expressed as 



1/2 -i^.0„/2 -i£TA„/2 



In terms of the new variables, Q thus reads 



sm 



[agfi/2] 



_-t/l/2giM03/2 

sm [cr(A^ - A„ +5M)/2] 



(117) 



(118) 
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Comparing (|93p and (IllSp one sees that for model RS the dual matrix matrix Q is obtained from L by changing 
a ■(r^ fj,, g ^ —g, p — ?> and q ^> A. It means that this model is self-dual (the matrices L and Q are the same up to 
a change of notation). One can show that the transformation from {qk,Pk) to action and angle variables {Xa,(f>a) is 
canonical ^2d\. As mentioned, a consequence of the unitarity of the matrix Q is that Va and Wa have the same sign 
as sin(A^T/2)/ sin(T/2) for all a. This implies that certain inequalities have to be verified by the A^, which we now 
derive in a way similar as in the previous section. 
Let us define the function 

h{x) = ^y^cot [{x - crX^)/2] ~tcot{NT/2) , (119) 

7 

which is periodic with period 27r and can be considered as a function on the unit circle. It has N poles a,t x — aX^. 
Taking the imaginary part of Eq. (|107|) . with c given by (jll2|) and given by (|117|) . one obtains that h{x) has N 
zeros at a; = aXj — t. When all V-y are positive, the same arguments as in the previous section imply that between 
two consecutive poles of h{x) there must be exactly one zero. It means that between two nearby eigenvalues aXa 
there is one and only one number of the form aX^ — r. A similar reasoning starting from matrix R = leads to 
the conclusion that between two consecutive aX^ there must also be one and only one number of the form aXa + r. 
These two conditions (shift by +r or — r) are equivalent, thus we can restrict ourselves to a shift by +t, i.e. the 
condition that the sets {cA-y, 1 < 7 < N} and {aX~^ + t, 1 < 7 < N} intertwine on the unit circle. It is a necessary 
and sufficient condition for the matrix (I102p to be an unitary matrix. A similar conclusion is readily obtained when 
all Vj are negative. 

The conditions implied by this type of intertwining have been discussed in [l^, HI]. There is an fundamental 
difference between these eigenvalue conditions in the RS model and in the CMt model discussed in the previous 
section. In the CMi case the Lax matrix is Hermitian, and the Aq. can take values on the whole real axis. In the RS 
case, as the Lax matrix is unitary, the aXa lie on the unit circle. Therefore poles and zeros cannot be ordered in a 
simple way as in the previous case, and the analysis of the previous section does not apply. 

For completeness we shortly repeat the arguments of the papers [2?, '28^ . Let us put all eigenvalues crXa of the Lax 
matrix (j93p on the unit circle and divide the circle into sectors with angle t. Denote the (positive) angular distance 
from the boundaries of the kth sector in counter clockwise direction by Xk and in clockwise direction by yk, as in 
Fig. [S] After a shift by t, the intertwining relations imply that only one of the two points corresponding to Xk and yk 







/ N't 




( V 



FIG. 5. Division of the circle into sectors of angle r. Only sectors numbers k, k + 1, and k + 2 are indicated. Black circles 
indicate the positions of the aXa, where exp(iaAa) are the eigenvalues of the RS Lax matrix. 

will fall in-between points corresponding to Xk+i and yk+i- The first case corresponds to Xk > Xk+i and yk+i > yk- 
In the second case the inequalities are reversed and Xk < Xk+i and yk+i < Vk- In both cases the inequality 

iUk+i - yk){xk+i - Xk) <0 (120) 

is fulfilled. 

Let us consider consecutive sectors of angle t as in Fig. [5l Denote the number of eigenvalues in each sector by 
Hk- After a shift by r, eigenphases from the fcth sector will move into the (fc + l)th sector. The Uk shifted points 
divide this sector into Uk + 1 intervals. As was proved above, eigenphases in the {k + l)th sector have to intertwine 
with these shifted eigenphases. Therefore all -I- 1 intervals except the first and the last will be occupied. The first 
will be occupied if Xk+i < Xk and the last interval will be occupied provided yk+i > yk+2- These statements can be 
rewritten in the form of the recurrence relation 



Uk+i =nk-~l + Q{xk - Xk+i) + Q{yk+i - yk+2) , 



(121) 
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where 8(t) is the Heaviside step function, Q{t) = 1 when t > and 8(t) = for t < 0. From (I120|) it follows that 
this relation can be rewritten in the form 

nk+i = n.k - 1 + Q[xk - Xk+i) + 'd{xk+2 - ifc+i) • (122) 
We now specialise to the case where r depends on the size N of the Lax matrix. We set 

27r 

r=— a (123) 

with fixed a. The case of integer a is trivial and we thus assume that a is not an integer. The total number of sectors 
of angle r in the unit circle is 



K = 



(124) 



where [t] denotes the integer part of t. Suppose the beginning of the first sector lies at position crXi. We choose to 
consider that this eigenvalue does not belong to the first sector, i.e. yi = 0. Applying Eq. (|120p for fc = implies that 
xi > xq. Thus necessarily n2 = ni + 1 and for A: > 3 one easily gets from Eq. (|122p that 

Uk = ni + Q{xk+i - Xk) ■ (125) 

The total number of eigenvalues lying into all K sectors obeys the inequalities 

K 

N -ni ~1 <^nk < N -1 . (126) 
fe=i 

The right-hand side inequality comes from the fact that when a is not an integer the union of all K intervals does not 
overlap the whole circle: in particular, it does not contain the first eigenvalue aXi from which we start our sectors. 
The left-hand side inequality is a consequence of the fact that if eigenvalues were shifted by — r, the first sector in the 
opposite direction would have exactly the same number of eigenvalues as the second sector, i.e. ni + 1 eigenvalues: as 
all K sectors does not cover the whole circle and the uncovering region is smaller than the sector of angle r, it follows 
that the number of eigenvalues in all K sectors is larger than N — (ni -|- 1). 
From (|125p one easily obtains a second inequality 

K 



Kni + 1 <^nk < K{ni + 1) . (127) 

From these two inequalities it follows that 



^ l<n,<^. (128) 



By definition of K we have 



K + 1 - - - X 



N N 

KK < — . (129) 

a a 



Substituting in the right-hand side of ()128|) the minimum of K and in the left-hand side the maximum value of K 
one gets 

a(a — 2) , , 

a-l-——<n,<a + ^ 130 

N + a N — a 

which entails that for N large enough a — 1 < ni < a. Since ni has to be an integer, ni = [a]. This means that for 
sufhciently large TV, within an interval of length r from any eigenvalue there are always exactly [a] eigenvalues. 

The inverse statement is also true. If within an interval of length r from any eigenvalue there exist exactly [a] other 
eigenvalues then all Va and Wa have the sign of sin(T/2)/ sin(iVr/2). To see this let us consider e.g. Va given by 
Eq. piop . It is a product of terms sin x/ sm{x — r) where x = a{Xa — A^)/2. As the distance between two eigenvalues 
may be restricted, < aXa — crXp < 27r, x obeys inequality < a; < tt. Therefore sin(a;) > 0, and sin(a; — r) is 
negative when < x < t and positive when t < x < tt. If within an interval of length r from aXa there are exactly [a] 
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eigenvalues, then in the product formula for Va there are exactly [a] negative terms, so that its total sign is (— 1)[°1. 
For N large enough the sign of sin(T/2)/ sin(A^r/2) with r = 2Tra/N is the sign of siuTra, which is precisely (— 1)[°1. 

The above arguments prove that for sufficiently large N (whose value depends only on a) the necessary and sufficient 
condition for the unitarity of matrix Q is that at distance |t| from any eigenvalue there exist [a] other eigenvalues. 

As matrices L and Q are dual, the unitarity condition for L can be readily deduced: it is that at distance |r| 
from one coordinate iiqk there are exactly [a] other coordinates. Note that in [20] only the case < a < 1 had been 
considered. 

These restrictions determine the allowed region in coordinate space. We choose as the 'natural' measure of momenta 
and coordinates the uniform distribution for momenta (between and 27r/cr) and coordinates uniformly distributed in 
the allowed region as explained above. After the change of variables from coordinates and momentum to action-angle 
variables it follows that the resulting distribution of eigenvalues will be also uniform but in the allowed region of 
eigenvalues with the only restriction that any interval ](TXa,aXa + 2-Ka/N[ contains exactly [a] eigenvalues. In next 
section we will show that it is possible to calculate asymptotic expressions for the joint distribution of eigenvalues by 
using a transfer operator technique. 

For numerical investigations, we consider an ensemble of unitary matrices of the form (j97|) . with pk chosen as 
independent random variables uniformly distributed between and 27r and the picket fence distribution of coordinates 
qj = j, I < j < N (see section |T|. Choosing constants such that /i = 27r/iV, a — 1, and g — a, so that r — 2Tra/N, 
direct calculations yield 



Vk = Wk 



siniVr/2 



iVsinT/2 

With a slightly different choice of phases in ([M]) , matrix L simplifies to 



(131) 



" ~ 1 - e2i^((fe-'-+'i)/^') ' 

where we denote pk = ^k- This is a particular specialisation of the Lax matrix for the Ruijsenaars-Schneider model. 
In the form (|132p but with a = bN with fixed b it first appeared in [23] as a result of the quantisation of a classical 
parabolic map on the torus proposed in (lo| . When 6 is a rational number the map considered in [lo| corresponds to 
a pseudo-integrable map of exchange of two intervals. For this particular case where a depends on N, the spectral 
statistics of the unitary matrix (|132p has been obtained analytically in [l^j and [2^ without knowledge of the relation 
with the Lax matrix of the Ruijsenaars-Schneider model. In the present case, where a is a fixed parameter independent 
on N, results are quite different. Analytical calculations of spectral correlation functions for this model are performed 
in the next sections. 

To illustrate the accuracy of the analytical results that we derive in the next sections, we show in Fig. [6] results of 
numerical computations for matrices of the form (jl32p for different values of the parameter a. Agreement is remarkable 
for all parameter values. 



VI. JOINT DISTRIBUTION OF EIGENVALUE SPACINGS FOR MODEL RS 

We now calculate asymptotic expressions for the joint distribution of eigenvalue spacings for the Lax matrix ensemble 
corresponding to the RS model. As discussed in the previous section, eigenvalues are such that within an interval 
of length a from any eigenvalue there are exactly [a] other eigenvalues. We introduce the rescaled nearest-neighbour 
spacings 

N 

4 = — (A„+i - A„) . (133) 

ZTT 

The answer strongly depends on the integer part of a. Therefore we consider different cases separately. 



A. < a < 1 



The simplest case corresponds to < a < 1. In this case the only restriction is that the distance between the 
nearest eigenvalues is larger than a, namely ^k > a. It is convenient to rewrite this restriction as follows. The joint 
probability density of having -I- 1 eigenvalues inside an interval of length L is given by 



^(6, 6, . . . , C^) = -T^TT-^ TT 9i^i)S [L-Tjk] (134) 
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FIG. 6. Nearest-neighbour distributions P{n,s) for the random matrices of model RS with fj, — 2n/N and g = 1/2 (top left), 
6/5 (top right), 4/3 (bottom left), and 9/4 (bottom right), averaged over 1000 realisations of matrices of size TV = 701. Solid 
lines are numerical results, dashed lines indicate the analytical expressions (|170|l for g = 1/2, (|176p - (ll80p for g = 1.2 and 
g = 4/3, and (PM)) - (fTMl) for g = 9/4. In each panel (from left to right) P(s) = P(l, a) (red), P(2, s) (green) and P(3, s) (blue). 



where 



and Zm{L) is the normalisation constant 



. ^ _ 1 when x < a 
9\^i \ I otherwise 



-^0 j=l V fc=i / 

We are interested in the joint probability distribution of n consecutive spacings 

/"OO /"OO 

P(Cl,---,fn) = / d^«+l---/ dCArp(Cl:6, ■ • ■ ,Cjv) 

Jo Jo 



(135) 



(136) 



(137) 



when L, — > cxD with mean level spacing A — L/N remaining constant. In the following we set A = 1. We shall 
proceed as it was done in [2^. The multiple integrals in (I137|) are easily calculated by introducing the function 



hn,N{L) 



3 = 1 \ k=l ) 



(138) 



whose Laplace transform reads 



5n.A^W = A(<)^-"n5(aK*^'= 



(139) 



fe=i 



Here 



A(<) = / 5(a:)e^'^da; : 



(140) 
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is the Laplace transform of g{x). The inverse Laplace transform of X{t)^ " is then 

-1 /•C+'ioO -1 /•C+'icxD -1 

— / A(0^-" exp(LO dt = — TTTt;: (^(In A(t) + At)) dt. 



(141) 



The large- A'^ behaviour of p4ip is obtained by saddle-point approximation. The value c corresponds to the solution 
of the saddle-point equation 



(recall that A = 1). The solution reads 



l + ^^O (142) 



c=-^. (143) 
1 — a 



Similarly, the inverse Laplace transform of gn.N{t) can be calculated by saddle-point approximation. Rather than 
calculating explicitly all prefactors coming from the integration, it is easier to observe that the large- behaviour of 
the normalisation factor = ho^N{L) is obtained similarly. One finally gets 

na/{l~a) " 

B. 1 < a < 2 

We now consider the case 1 < a < 2. According to previous section for sufficiently large matrix size, any eigenvalue 
Xk is such that there exists exactly 1 eigenvalue in the interval ]xk,Xk + 2'Ka/N[. In other words, the constraints on 
Xk are that Xk+i &]xk,Xk -t- 2Tra/N[ and Xk+2 > Xk + ina/N. In terms of the differences (|133p between consecutive 
eigenvalues, the above restriction are equivalent to 

< Cfc+i < a , a < Cfc+i + 6 • (145) 

Introduce the function f{x) by 

, , _ J 1 when < X < a ^^^(\\ 

•' ^ ' ~ [ otherwise ^ ^ 

and g{x) as in (|135|) . Then the joint probability density of A^ + 1 eigenvalues inside an interval of length L is given 
by the following expression 



jv / ^ \ 



p(?i, 6, .■.,^n) = 1 1 fi^M^j + G+i)<5 [L-yjk] (147) 

where Z]\j{L) is the normalisation constant 

Zn{L) = / d^i . . . / dCw I I /(0)5(e, + ^j+i)S [L-yjk] ■ (148) 



poo poo N / N \ 

/ . . . / dCw n /(o)5(e, + ^j+i)s u - E 

•^0 -^0 j=l \ k=l / 



The large- A^ behaviour of the joint probability distribution of n consecutive spacings (jl37p can be then obtained as 
above, following [29| . We review the main steps of the procedure. Introducing the function 



poo poo N I ^ \ 

\ de„+i . . . / d^AT [] /(e,).g(G + o+i) M ^ - E ' 

•^0 -^0 j = l V fc=l / 



hn^N{L) = / de„+i . . . / d^AT I I /(e,).g(G + O+i) M ^ - > . ?H ' (149) 
its Laplace transform reads 

N 

9n.N(t) - / de„+i . . . / dCjv TT e"'*V(e,)5(e, + O+i) • (150) 



nOO r-OO 

/ de„+i . . . / dCjv n e"*«V(e,)5(e, + O+i) 
Jo Jo j^-^ 



22 



This quantity can be seen as a product of transfer operators 

ff„,jvW = -?ft(a,6)^t(6,6).-.i^t(?n-i,60 (151) 

X de„+l.../ diNKt{^n,^n+l)--.Kt{^N-uU}Kt{iN,^l), (152) 

Ja Jo 
where the transfer operator, Kt{^,^') is defined as 

Ktit n = /(05(e + ana e-*(«+«')/2 . (153) 

This is a real symmetric operator. Its real eigenvalues, Aj(i), and eigenfunctions, (j)j(t;£^), verify 

Kt{^,a<t>At-^ad^' = A,(i)0,(<;C) • (154) 
As this operator is real symmetric its eigenfunctions can be chosen to be orthonormal 

^,it;OMt;Od(^S,k . (155) 
The transfer operator can be expanded over the basis of eigenfunctions as 

3 

In the large- iV limit the dominant contribution to (jl51l) is given by the largest eigenvalue Xo{t). Then using the 
orthogonality property (|155p of the we get 











9n,Nit) ~ Y[Kt{^k,^k+i)Mt)''-''^'Mt;UMt;ii) ■ (i57) 

fc=l 

The large-TV behaviour of /in,jv(-^) is obtained by performing the inverse Laplace transform of \o{t)^~^'^^. As in 
(|141l) the leading term is obtained by saddle-point approximation; here the saddle-point is such that 

An(c) , , 

Ao(c) 

(again we take A = 1). Calculating in a similar way the large- A'' behaviour of the normalisation factor Zn = /io,Ar(i), 
one finally gets 

p(ei,6,---,6O = T;r^'/'o(c;a)ifc(ei,e2)i^c(6,C3)...i^c(e„-i,e«)0o(c;e«) . (159) 
\ (c) 

C. m < a < m + 1 

The general case m < a < m+l can be treated in a similar way. In this case each interval ]xk,Xk + 27ra/N[ contains 
exactly m eigenvalues Xk+i, ■ ■ ■ , Xk+m- In terms of the rescaled differences (|133p between consecutive eigenvalues this 
condition is equivalent to the following two inequalities 

< 6 + Cfc+i + • • • + a+m-i < a (160) 

a < Cfe+l + • • . + ik+rn ■ (161) 

In analogy with Eq. (jl47p . the joint probability of eigenvalues spacings thus reads 

TV / ^ \ 

p(Cl, 6, • • • , Cjv) - n /(^^- + ■ • ■ + 0+™-l)5(0 + • ■ ■ + ^J+rn)S U - E (162) 

' 3=1 \ k=l I 
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with / and g defined by (I146P and (|135p . The large- iV behaviour is calculated as above by introducing a transfer 
operator, which in this case depends on two sets of variables | = (Ci, • • ■ I'Cm) and ^' = (CJ, ■ ■ ■ ,^^) shifted by one 
unit i.e. ^2 = Ci, = ■ • Cm = Cm-i (see e.g. [29]). The explicit form of the transfer operator is the following 

m.i')^5{^2-^[)...5{u~i'm-i) 

xe-*«i/V(6 + • ■ • + U)9{£.i + + OfiCi + ■■■ + Oe-*^"/' • (163) 
The eigenvalue equation 

'mO</'(€')dr = A0(O (164) 



/- 



reduces to a one-dimensional equation because of the ^-functions appearing in the definition of the transfer operator. 
This equation can be written in the form 

POO 

g-t6/2 / e-''/^gi^,+C2 + ---+Cm + zm;^2,.--,U,z)dz = Ximt;^,,...,U) ■ (165) 
Jo 

Here it is implicitly assumed that all variables > and 

(f){t; Ci, . . . , Cm) = when £.1 + . . . + U > a . (166) 

As above the largest eigenvalue, Xo(t) as well as the corresponding eigenfunction (f>o{t;^), calculated at point t = c 
obeying the same saddle-point condition Eq. (|158l) . determine all correlation functions in the limit of large N. The 
joint probability of n consecutive spacings takes a different form for n < m and n > m. For n < m 



j*a pa 
P(Cl,---,Cn) = / ••• / 

JO Jo 



d6n0o(c;6, ■ • ■,£,ni)(t>o{c]£,m,-- . ,fl), (167) 



while for n > 



. . . ,60 = Ao(c)-"+™^o(c;e„,. . . ,6.-™)</>o(c;6, • • • ,U)e-^^"-«= 

n— m+1 n—rn 

X n /(Cj- + • • • + cj+™-i) n 3(0 + ■ • • + (168) 



VII. NEAREST-NEIGHBOUR SPACING DISTRIBUTIONS FOR MODEL RS 

In the previous section we have derived expressions for the joint distribution of eigenvalue spacings p{S,i, 
From these expressions the nth nearest-neighbour spacing distribution can be calculated as 



(169) 



A. < a < 1 

For < a < 1 these integrals are easily calculable (e.g. by Laplace transform) and from the joint distribution 
Eq. (|144p we obtain 

( e("''-'*)/(i-'*) (s - na)"-i 
P{n,s)^} (l_a)" (n-1)! ' « > . (170) 
[ 0, < s < na 

Comparison with numerical simulations is displayed at Fig. [51 
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B. 1 < a < 2 



When 1 < a < 2 the joint distribution is given by Eq. (|159p . What remains is to calculate the largest eigenvalue of 
the transfer operator i^t, as well as its associated eigenfunction. As Kt is a positive operator, the analog of the Perron- 
Frobenius theorem states that the eigenvector corresponding to the largest eigenvalue is positive. Orthogonality of the 
eigenfunctions implies that the converse is also true. Thus if one finds a positive eigenfunction then the corresponding 
eigenvalue is the largest one. The eigenvalue equation (|154p is equivalent to 

e-*«/2 r e-'i'/'mdC = Xm ■ (171) 

Let us look for solutions of Eq. (I17ip positive on [0, a] under the form (/)(^) = sinhpf , with p some unknown complex 
parameter. Since Eq. (|17ip should hold for all f G [0,a] we get the necessary condition 

t = -2pcoth{pa). (172) 

When t < —2/a, Eq. (|172p admits two real solutions p ~ ±pa. Thus (/){£,) — sinhpo'f with po > solution of Eq. (|172p 
is a positive solution of Eq. (|17ip . li t > —2/a, Eq. (|172p admits two pure imaginary solutions p = ±ipo, thus 
4>{^) = sinpoC with po > and ipQ solution of Eq. (|172p is a positive solution of Eq. (I17ip . Finally ii t — —2/a, p = 
is the unique solution to Eq. ()172p . In that case (/)(^) = ^ is a solution of Eq. ()17ip which is positive on [0,a]. Thus 
for all t we have a positive solution to Eq. (|17ip . Properly normalised, this solution gives the eigenvector (j)o{t;^). 
The corresponding eigenvalue is given by 

„(p-t/2)Q 

Mt) = T7:r, (173) 

p-t/2 

with p an implicit function of t. 

The saddle-point c is a solution of Eq. (|158p . For Aq given by Eq. ()173p the condition becomes 

1 + 2a(2 - a)p^ - cosh 2pa + 2p{a ~ 1) sinh 2pa ^ (174) 

and the saddle-point c is obtained from p through Eq. (I172p . Equivalently, this condition can be expressed as 

22:^-zsinh2z . 

a = ; — 2 : J z = pa . (175) 

-\- sinh z — z sinh 2z 



Q 1.5 



ap 



FIG. 7. (Color online). Graph of Eq. (|175p for real (dashed) or pure imaginary (solid) z — pa. Dashed vertical line indicates 
the abscissa equal tt. 

In Fig. [7] we plot a as a function of z = pa. For 1 < a < 4/3 Eq. ()175p has a unique real solution po > 0, 
and 00(0 = sinhpoC is a positive eigenfunction of the transfer operator. For 4/3 < a < 2 Eq. (|175p has a unique 
pure imaginary solution ipQ with po > 0. Furthermore in that latter case poa G]0,7r[, so that <poiO — sinpo'? is 
an eigenfunction of the transfer operator which is positive on [0,a]. At a = 4/3 the unique solution is p — and 
= C is a positive eigenfunction of the transfer operator. The nth nearest- neighbour spacing distribution can 
now be calculated from Eq. (|159p . 
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In the case n = 1 it directly gives us the nearest-neighbour spacing distribution P{s) = A^0o(s)^, where A is the 
normahsation constant. It is nonzero only for s £ [0, a], where it takes the following form 

{A^ sinh^(/9s) when 1 < a < 4/3 
fis2 when a = 4/3 . (176) 

sin^(/9s) when 4/3 < a < 2 

Constants A and p can be determined either by solving Eq. (|175p and normalizing the eigenfunction 0o(C)i or equiv- 
alently by imposing the normalisation conditions The next-to-nearest distribution, P(2, s) is non-zero only when 
a < s < 2a and within this interval it is given by 

P(2, s) = T^^e-^^/^ r MOM-^ - m ■ (177) 



In particular for a = 4/3 all integrals can be calculated analytically and P{2, s) has the form 

P(2,.) = (-^ + |.-^.V^/4-i. (178) 

In a similar manner one can obtain the higher nearest-neighbour functions. For example, ^(3, s) is given by the 
formula 

a2 -cs/2 j-a /.a , /"^ / " \ 

ms) = dei</)o(6) j '^^^'^ J d60o(S3)<5 ( ^-E^O • (179) 

It is non-zero only when a < s < 3a. In particular for a = 4/3 we obtain 

v(3s)-f (t-i^+^^Vf^-' + i^; when 4/3 <.< 8/3 

^_9^27^,_^^3)g3./4-l_^gg3./2-4 when 8/3 < S < 4 " ^^^"^ 

In Fig. [6] these formulas are compared with numerical simulations and show a remarkable agreement. 

C. 2 < a < 3 

For 2 < a < 3 the joint distribution is given by (I168p . The largest eigenvalue and corresponding eigenfunction of 
the transfer operator (|163p arc solution of the eigenvalue equation (|165p . For m = 2 it takes the form 

e-*«^/' r e-*«^/2^(e2,6)d6 = A0(a,6) • (181) 
Let us look for solutions of the form similar to Bethe Ansatz 

<^(^l,^2) = e"^i+^^" ^Q-Kl + ia~P-I^K2 _|_ g(-a-(-/^-(-M)?l-a?2 

where we have set fj, = —t/2. As Eq. (|18ip has to be fulfilled for all ^1,^2, this function is a solution of (I18ip if and 
only if the following conditions are valid 

pa(/j+/3) pa(p-a) „a{a-l3) \ 

p, + p fi ~ a a~p a 

From the first equality in Eq. (|183p one can express /i as a function of a and /3. After inspection we found that the 
solutions of the above equations have the following form 

1^1 1 , 

aa = —xi + 1x2 , pa = "2^-*^ ' ^'^ ~ ~^ (184) 

with real parameters Xi, X2, and X3. Under this substitution the eigenfunction ()182p is transformed to 



<A(ei,6) = e-^«^-«^)/2« (^sin^^^^^±^ -e(— ^^^^/-^sin^ -e-(— ^)«^/"sin^ 



(185) 
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where from (I183|) Xi, X2, and 2:3 must fulfill the following equalities 

xi X3 + 1x2 X3 - 1x2 a 
and depend on time t through the relation 



(186) 



This implies that 



ta Xi , 
y = y + -3. (187) 

LcLll X2 

and, consequently, xi is related with X2 as follows 

— = !i!i^e"^/*^'^"^ . (189) 

Xl X2 

The eigenfunction corresponding to the largest eigenvalue of the transfer operator is thus given by (|185p with xi, X2, 
and X3 real parameters depending on t, which must verify (|188p and (|189p and be such that (f>{£,i,^2) is a positive 
function over [0, a]^. 

The saddle-point condition is again given by Eq. (|158p . Using (|186p - (|189p we get a second relation between xi and 
X2, namely 

1 , 2-sin(2a;2)/a;2 ^^^^^ 



1-1/xi 1 + sin^(x2)/a;2 - sin(2a;2)/a;2 

Equations (I189P and (I190p determine parameters xi and X2 at a given a. To get a positive eigenfunction (f> it is 
necessary to get the solutions in the intervals 

Xl <0, TT < X2 <27r . (191) 

The knowledge of these parameters allows us to calculate the eigenfunction (|185p , from which the nearest-neighbour 
distributions can be deduced through Eq. (|168p . The first distributions read 

P{s)^A 0(s,2/)0(y,s)dy , (192) 
Jo 

P{2,s)^a( ^{s-y,y)(t>iy,s-y)dy , (193) 
"'0 



and 



A 



S — X 



P{3, s) = T- / dxe^'' / d?/e^2/0(s ~^x~y, x)4){s ~x~y, y)dy , (194) 

^ J s — a J s—a 

with A the normalisation constant 

A=(^j^dx(^j^ Ay^{x,y)<j,{y,x)y^ . (195) 
These analytical expressions perfectly agree with numerical simulations, as shown in Fig. [51 

VIII. LEVEL COMPRESSIBILITY FOR MODEL RS 

The expressions for the joint distribution of eigenvalue spacings p(^i, . . . , fn) obtained in section IVTl allow to derive 
formulas for the level compressibility Xj which characterises the asymptotic behaviour of the number variance. 

The number variance Y?{L) is the average variance of the number of energy levels in an interval of length L. It is 
defined from the two-point correlation function R2{s) = X]n=i ^("•' 

I]2(L)=i-2/" ds{L~ s){l- R2{s)). (196) 
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For systems with intermediate spectral statistics, ~ xL for large L. In order to obtain the large- iV behaviour 

of the level compressibility we calculate the Laplace transform of the two-point correlation function. It has a series 
expansion of the form 

32 W = ds R2{s)e-'^ = j + ^ + 0(t) (197) 
which allows us to obtain x (see [29] for more detail). 



A. Case < a < 1 



The nth nearest-neighbour spacing distributions are given by Eq. (|170l) . Summation over n gives the two-point 
correlation function, and its Laplace transform is readily obtained, yielding 

- e-^(l+/-at)-l - ^'''^ 

Small-t expansion of 52 (t) gives 

X = (!-«)'■ (199) 



B. Case 1 < a < 2 



= ^n-l, ^ I d6 • • ■ / d^„ 



The functions P(n, s) are given by Eqs. (|159p and (|169p . Their Laplace transform reads 

1 

Ao~(c) ^0 

X </.o(c;6)ifc(6,6)A'c(6,6)..-i^c(e.-i,C«)'/'o(c;en)e-*(«^+-+«"), (200) 

where as in the previous sections Aq (c) is the largest eigenvalue of the transfer operator (|153p and (f)^ (c; ^) its associated 
eigenfunction, both taken at the saddle-point c. Using the definition (|153p of the transfer operator, we see that g{n, t) 
can be rewritten 

9M = TK^J / de'0o(c;Oe-*«/'ifc+t(?,n""Vo(c;e')e-*«'/^ (201) 
Ag (c) Jo Jo 

Replacing the transfer operator by its expansion ()156p and summing over n we get 

One can check, using normalisation (|155p of the eigenfunctions and the saddle-point condition (|158p . that the leading- 
order term is given by 32 (0 ~ 1/i- The next-order term can be simplified using the normalisation of It yields 



from which one gets 



Here Ao(i) is given by (|173p (with p depending on t through (11721) ). and c is given by condition (|158p . After calculation, 
X can be expressed as a function of p at the saddle-point. We get 

fa? 4a(l — a)z^ + sinh^ z ,9 \ sinh^ z 
X = ,^ ■ , o ^9 smh^ z z = pa, (205) 



4 (2z-sinh2z)2 

with p the real positive solution po of p75p for 1 < a < 4/3 or the pure imaginary solution ipo of (|175p for 4/3 < a < 2. 
For a = 4/3, the limit p ^- in pI5)l gives x = 4/9. 
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C. Case 2 < a < 3 

As in the previous case x is given by (|204p with Ao(t) given by p86p . with xi, X2, and t related through 
(|187p - (|189p . From (|187p - (|188l) . parameter xx can be expressed as 

XI = -2— at. (206) 

tan X2 

Differentiating both (|189p and (|206p with respect to time we obtain Axi/dt and Ax2lAt as a function of xi and X2, 
and then similarly d^xi/dP and d^a;2/dt^. Using (|186l) . the saddle-point condition (I158P can be rewritten 

and from (|204l) x can then be expressed as 

Using the expression obtained d^xi/dt^ we finally obtain x as a function of xi and X2, with xi, X2 obtained as solution 
of (Hinil-dTnni). inverting we get 



3in^ X2 + {a — 2)x2 + (1 — a)x2 sin 2x2 
(a — 1) sin^ X2 + {a — 3)x2 + (2 — a)x2 sin 2x2 



_ Lioiii .1.2 ~r (u, — ^ju,2 -r y± — u,)a,2 aiii-^.1,2 ^'onn^ 

— 7 ^ , . 2 } , 2 . /.^ ^ '■ — ^ — ■ (,^uyj 



After some manipulation x simplifies to 

X = , . 2 _^ 2 r^--^[(a-3f(a-2)x^ (210) 

a(sm X2 + ^2 — X2 sm 2x2) 

- (a - 3)(a - l)(2a - 5)x^ sin2x2 + 2(a - 2)((cos2x2 + 2)(a - l)(a- 2) - 3)x^ sin^ X2 

— 2a{a — 2) (2a — 3)x2 COSX2 sin'^ X2 + a{a — 1)^ sin^ X2] . 

Figure |8] is a plot of the level compressibility x- The theoretical prediction obtained from (|199p . (|205p and (|210p 
agrees with numerical data. 

D. Asymptotics in the vicinity of integer a 

For integer a the spectrum is rigid and thus the level compressibility is expected to take the value 0. Here we 
consider the first-order expansion of x in tiis vicinity of integer a. We will show that at lowest order the expansion of 
X around a = n is given by x — (1 ~ a)^/n^. 

We first consider the expansion around a = 1. Let a = 1 + e. For a < 1 we have X = (1 — a)^ = e^, thus expansion 
is trivial. For a > 1 x is given by (|205p with a and z related by (|175p . At a = 1 the solution of (I175P is 2; = 00. An 
asymptotic expansion of pOSp and (|175p yields 

+ F{z)e-^' + o{e-^') (211) 



2z - 1 

where F{z) is some rational fraction in z. Thus at first order 



2 V e 

Expanding x for large z to lowest order gives 



z=Ul + -]. (212) 



X = a{a-l) + a^ (^-^] + G{z)e-^' + o{e-^') (213) 



.4z^ 2z, 

with G{z) is some rational fraction in z. Using (|212l) one gets x — = (1 — o)^ 
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FIG. 8. Level compressibility x- Black diamonds are the numerical values extracted from a cubic fit of over the range 

L € [0,80], with T?(L) the variance of the number of levels in an interval of length L averaged over 50 windows of length L, 
calculated from the unfolded spectrum with mean level spacing A = 1, for 10000 realisations of the random matrix with matrix 
size iV = 256. Red curve is the theoretical prediction (fTWl) . pOS]! and (pTOl) . 



Suppose now that a = 2 + e. For a < 2 x is given by (j205p . with p = ip^ solution of (|175p . Equivalently, x is given 

by 

/ Aail — a)z'^ + a? svc? z , \ sin^ z 
X= ^ + ^77Tr^-3T7T2 sin^ (214) 



(215) 



4 (2z-sin2z)2 
with z the real positive solution of 

2z^ — z sin 2z 



— 2z sin 2z + sin^ z 



At point a = 2 the solution is z = tt. Expanding both sides of (|215l) at lowest order in e with z = tt + zie we get 
zl = 7r/2. Inserting this expansion for z in (|214p we obtain x — £'^/4 + o{e^). 

For a > 2 X is given by (I210p . with xi and 2:2 specified by (|189l) - (jl90p . At a = 2, we have 



2 sin 2:2 — X2 sin 2a;2 
2(2:2 + sin^ X2 — X2 sin 2x2) 

which vanishes for X2 — tt. For X2 = n + tie + t2€^ we have the expansion 



X = 2 I -2 r^;-^ (216) 

" ---^ _L c'Tvi'^ -- ' -- ' 



X = - ^) ^ + f4 + 1^ - ^ - + «(^^)- (217) 



Equation ()189p is equivalent to 



X2 \ xisina;2 . . 

exp xi - = (218) 

imiX2 J X2 



and the small-e expansion of both members of this equation reads 

xi-^ = --^ + ^-l + o(l), (219) 

^2 ^ .^t,-5.tf + 6t3 ^^t2-4.t,^2^3 + ^(^4) (220) 



2:2 _ 


TT 


7ri2 


tan 2:2 






xi sin a;2 




2i2 


0:2 


7r2 
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(we use (|209p to obtain the expansion of ii). This iinphes that the two first terms in the expansion (|220p must vanish, 
thus ti = n/2 and t2 = 0. Putting these values into (|217l) gives x = + o{e^). 

The same result can be obtained from (|210|) and (|189|) ~ (|190|) for a < 3. At a = 3 again X2 takes the value tt, and 
an expansion X2 ^ tt + tie gives ti = tt/3, whence x = + o{e^). 



IX. CONCLUSION 



In this paper we construct new random matrix ensembles with unusual properties. Random matrices from these 
ensembles are Lax matrices of TV-body integrable classical systems with a certain measure of momenta and coordinates. 
Though such matrices are not invariant over rotation of the basis (as usual random matrix ensembles) the joint 
distribution of their eigenvalues can be calculated analytically. Four different models are considered in detail. Three 
of them correspond to rational, hyperbolic, and trigonometric Calogero-Moser models. The fourth is related to 
the trigonometric Ruijsenaars-Schneider model. For the trigonometric Calogero-Moser model and the Ruijsenaars- 
Schneider model spectral correlation functions are calculated explicitly. For rational and hyperbolic Calogero-Moser 
models Wigner-type surmises are proposed. Our formulas are in a good agreement with results of direct numerical 
calculations. 



Appendix A: Hamilton-Jacobi equations 



In this appendix we check that the action-angle variables Aq and (t>a for model CM,, verify Hamilton-Jacobi equations 
by calculating their time derivative. We use the fact that for the Lax pair (L, M) the matrix M can be seen as a 
time derivative operator for the eigenfunctions of the matrix L. Namely, if {uk)i<k<N is a normalised eigenvector of 
L, then 

uii. — Mf^r^r and = — U*Mrk {^^) 
r r 

(here * denotes complex conjugation). For model CM^, the Lax matrix M is given by 

Mkr = -igSkr^-, — 2 +i.g(l - 4,-)7 — - — ■ (^2) 
One can easily check that for 1 < fc, r < one has 

PrSkr + Mkriqk - Qr) ^ Lkr ■ (A3) 

Deriving p9l) with respect to time, using the definition of Q, yields 

(pa = ^[ul{a)qkUk{a) + ul{a)qkUk(a) + ul{a)qkUk{a)] . (A4) 
k 

From Hamilton-Jacobi equations qk = Pk- Using (jAip . we get that the time derivative of is given by 

= ^u*k{a)[pkSkr + Mkriqk ~ qr)]ur{a) ^^ul{a)LkrUr{a) = Xa ■ (A5) 

k,r k,T 

The time derivative of Aq is easily obtained from ([7l). ()Aip and (p4|) . yielding Aq = 0. This shows that the Aq and 4>a 
verify Hamilton-Jacobi equations. 



Appendix B: Identities 



The purpose of the Appendix is to give, for completeness, the proofs of certain often used formulas. 

Let coefficients hm obey the following system of linear equations for all n = 1, . . . , A^ with known Xm and ym 

m=l ^" 
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Then for all m = 1, . . . , N are expressed through Xm and ym by using e.g. the Cauchy determinants 



The following identities are also useful. For all / = 1, . . . , one has 



^ brn ^ UnMy^-Vn) , . 



N N 

^b„i=^{x,n-yni), (B4) 
m— 1 m— 1 

and 

m— 1 n 



A simple way to check (|B2p is to consider the function 

Asymptotically this function decreases as 1/x when a; — t- cxd, so that the integral over a large contour encircling all 
poles equals 1. Rewriting this integral as the sum over all finite poles gives 

^ _ ^ rirC^™ ~ yr) j-gy-j 

,^ i^m ~ Vn) Y\s^mi-^">- ~ -^s) 

which proves (IB2p . 

The equality (|B3|) can be obtained by the integration of the function 

over a contour which includes all poles. As this function decreases as 1/a;^ when x — ?> cxd the integral equals zero. 
Taking the sum over poles at x = Xn with all n = 1, . . . , A'^ and at x = yi one verifies (|B3|) . 
To get (|B4[) one has to integrate the function 

N 

/(^) = Yl (B9) 

n— 1 

over the large contour and compare the residues at infinity and eit x = x^. 
Let us now consider the function 

It decreases as 1/x when x oo and has poles at x = and x = Xm with rn = 1, . . . , A'^. Integrating it over a contour 
encircling all poles one obtains (jBSp . 
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